Predicting passsenger traffic EMT¶
Utilizaremos Prophet, una biblioteca para pronósticos de series de tiempo lanzada por Facebook y de código abierto el 23 de febrero de 2017. En este caso vamos a predecir el número predeciremos el numero de pasajeros de la linea 89
El pronóstico de series temporales encuentra una amplia aplicación en el análisis de datos. Estas son solo algunas de las predicciones concebibles de tendencias futuras que podrían ser útiles:
La cantidad de servidores que necesitará un servicio en línea el próximo año. La demanda de un producto de supermercado en un supermercado en un día determinado. El precio de cierre de mañana de un activo financiero negociable. Para otro ejemplo, podemos hacer una predicción del rendimiento de algunos equipos y luego usarlo como línea de base: primero para establecer objetivos para el equipo y luego para medir el rendimiento real del equipo en relación con la línea de base.
Existen bastantes métodos diferentes para predecir tendencias futuras, por ejemplo, ARIMA , ARCH , modelos regresivos , redes neuronales .
Importamos Librerias¶
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import datetime
import seaborn as sb
import warnings
import re
import unicodedata
from unicodedata import normalize
from datetime import datetime
from fbprophet import Prophet
import logging
logging.getLogger().setLevel(logging.ERROR)
warnings.filterwarnings('ignore')
pd.pandas.set_option('display.max_columns', None)
plt.rcParams['figure.figsize'] = (20,12)
pd.set_option("display.float_format",lambda x:str(round(x,6)))
Cargamos datos¶
Nos centraremos en la Linea 89 Circular ronda tránsits
uno = pd.read_table("carga2019-05-15.txt", sep= ",",encoding='ISO-8859-1')
dos = pd.read_table("carga2019-05-16.txt", sep= ",",encoding='ISO-8859-1')
tres = pd.read_table("carga2019-05-17.txt", sep= ",",encoding='ISO-8859-1')
cuatro = pd.read_table("carga2019-05-18.txt", sep= ",",encoding='ISO-8859-1')
cinco = pd.read_table("carga2019-05-19.txt", sep= ",",encoding='ISO-8859-1')
seis = pd.read_table("carga2019-05-20.txt", sep= ",",encoding='ISO-8859-1')
siete = pd.read_table("carga2019-05-21.txt", sep= ",",encoding='ISO-8859-1')
ocho = pd.read_table("carga2019-05-22.txt", sep= ",",encoding='ISO-8859-1')
nueve = pd.read_table("carga2019-05-23.txt", sep= ",",encoding='ISO-8859-1')
diez = pd.read_table("carga2019-05-24.txt", sep= ",",encoding='ISO-8859-1')
once = pd.read_table("carga2019-05-25.txt", sep= ",",encoding='ISO-8859-1')
doce = pd.read_table("carga2019-05-26.txt", sep= ",",encoding='ISO-8859-1')
trece = pd.read_table("carga2019-05-27.txt", sep= ",",encoding='ISO-8859-1')
catorce = pd.read_table("carga2019-05-28.txt", sep= ",",encoding='ISO-8859-1')
quince = pd.read_table("carga2019-05-29.txt", sep= ",",encoding='ISO-8859-1')
dieciseis = pd.read_table("carga2019-05-30.txt", sep= ",",encoding='ISO-8859-1')
diecisiete = pd.read_table("carga2019-05-31.txt", sep= ",",encoding='ISO-8859-1')
dieciocho = pd.read_table("carga2019-06-01.txt", sep= ",",encoding='ISO-8859-1')
diecinueve = pd.read_table("carga2019-06-02.txt", sep= ",",encoding='ISO-8859-1')
veinte = pd.read_table("carga2019-06-03.txt", sep= ",",encoding='ISO-8859-1')
veinteuno = pd.read_table("carga2019-06-04.txt", sep= ",",encoding='ISO-8859-1')
veintedos = pd.read_table("carga2019-06-05.txt", sep= ",",encoding='ISO-8859-1')
veintetres = pd.read_table("carga2019-06-06.txt", sep= ",",encoding='ISO-8859-1')
veintecuatro = pd.read_table("carga2019-06-07.txt", sep= ",",encoding='ISO-8859-1')
veintecinco = pd.read_table("carga2019-06-08.txt", sep= ",",encoding='ISO-8859-1')
veinteseis = pd.read_table("carga2019-06-09.txt", sep= ",",encoding='ISO-8859-1')
veintesiete = pd.read_table("carga2019-06-10.txt", sep= ",",encoding='ISO-8859-1')
veinteocho = pd.read_table("carga2019-06-11.txt", sep= ",",encoding='ISO-8859-1')
veintenueve = pd.read_table("carga2019-06-12.txt", sep= ",",encoding='ISO-8859-1')
treinta = pd.read_table("carga2019-06-13.txt", sep= ",",encoding='ISO-8859-1')
treintauno = pd.read_table("carga2019-06-14.txt", sep= ",",encoding='ISO-8859-1')
Concatenamos datos verticalmente ya que van del día 1 al 31 Mayo
df = pd.concat([uno,dos,tres,cuatro,cinco,seis,siete,ocho,nueve,diez,once,doce,trece,
catorce,quince,dieciseis,diecisiete,dieciocho,diecinueve,veinte,
veinteuno,veintedos,veintetres,veintecuatro,veintecinco,
veinteseis,veintesiete,veinteocho,veintenueve,treinta,treintauno], axis=0)
df.sample(5)
df.shape
Preprocesado de datos¶
Filtramos por la linea 89 y renombramos variable
linea89 = df[df['NOMBRE_ITINERARIO']=='Circular ronda tránsits'].copy()
linea89.sample(5)
Nos centraremos en la variable FECHA Y VIAJEROS(Validadciones):
Con esta funcion lo que hacemos es forzar un 'guion' ya que hemos tenido muchos problemas con el formato datetime. El problema era que al hacer directamente datetime sobre columna FECHA, habian algunos datos que me invertia el orden, ponia el mes en medio y a veces al final en otros casos. Asi que pensé esta solución.
# Año-dia-mes
def fechas_formato(row, var):
if row[var] == '01/05/19':
return "2019-01-05"
if row[var] == '02/05/19':
return "2019-02-05"
if row[var] == '03/05/19':
return "2019-03-05"
if row[var] == '04/05/19':
return "2019-04-05"
if row[var] == '05/05/19':
return "2019-05-05"
if row[var] == '06/05/19':
return "2019-06-05"
if row[var] == '07/05/19':
return "2019-07-05"
if row[var] == '08/05/19':
return "2019-08-05"
if row[var] == '09/05/19':
return "2019-09-05"
if row[var] == '10/05/19':
return "2019-10-05"
if row[var] == '11/05/19':
return "2019-11-05"
if row[var] == '12/05/19':
return "2019-12-05"
if row[var] == '13/05/19':
return "2019-13-05"
if row[var] == '14/05/19':
return "2019-14-05"
if row[var] == '15/05/19':
return "2019-15-05"
if row[var] == '16/05/19':
return "2019-16-05"
if row[var] == '17/05/19':
return "2019-17-05"
if row[var] == '18/05/19':
return "2019-18-05"
if row[var] == '19/05/19':
return "2019-19-05"
if row[var] == '20/05/19':
return "2019-20-05"
if row[var] == '21/05/19':
return "2019-21-05"
if row[var] == '22/05/19':
return "2019-22-05"
if row[var] == '23/05/19':
return "2019-23-05"
if row[var] == '24/05/19':
return "2019-24-05"
if row[var] == '25/05/19':
return "2019-25-05"
if row[var] == '26/05/19':
return "2019-26-05"
if row[var] == '27/05/19':
return "2019-27-05"
if row[var] == '28/05/19':
return "2019-28-05"
if row[var] == '29/05/19':
return "2019-29-05"
if row[var] == '30/05/19':
return "2019-30-05"
if row[var] == '31/05/19':
return "2019-31-05"
linea89['DATE'] = linea89.apply(lambda row: fechas_formato(row, 'FECHA'), axis = 1)
linea89[['FECHA','DATE']].sample(5)
linea89.dtypes
Ponemos bien los nombres de las variables ya que vemos que estan desplazadas
linea89 = linea89.rename(columns = {' VIAJEROS': 'VIAJEROS',' TRANSBORDOS': 'TRANSBORDOS', ' RUTA':'RUTA', ' PARADA':'PARADA'})
linea89.dtypes
Nos quedamos con las dos variables que vamos a utilizar, y vuelvo a renombrar la variable:
df89 = linea89.copy()
df89 = df89[['DATE','VIAJEROS']]
df89.sample(4)
Hemos visto anteriormente que la variable DATE esta como object, hay que pasarla a datetime:
df89['DATE'] = pd.to_datetime(df89['DATE'],format='%Y-%d-%m')
df89.dtypes
# Finalmente se queda AÑO-MES-DIA
df89.sample(10)
Predecimos con Prophet y para ello preparamos el dataset como una serie temporal.¶
Cuantos validaciones hubieron cada día?¶
df89 = df89.groupby(['DATE']).sum().reset_index()
df89.sample(5)
# Prophet nos pide que trabajemos con estos nombres de variables
df89.columns = ['ds', 'y']
df89.tail(n=3)
df89.head()
En esta primera parte ajustamos con todos los datos (no hay división de prueba y entrenamiento.)
m = Prophet()
m.fit(df89)
Predecimos los siguientes 10 días
future = m.make_future_dataframe(periods=10,freq = 'D')
future.tail(10)
Predecimos valores Prophet pasando las fechas para las que queremos crear un pronóstico.
yhat representa la predicción, mientras que yhat_lower y yhat_upper representan el límite inferior y superior de la predicción, respectivamente.
forecast = m.predict(future)
forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()
from fbprophet.plot import plot_plotly
import plotly.offline as py
import plotly.graph_objs as go
py.init_notebook_mode()
#m.plot(forecast);
fig = plot_plotly(m, forecast)
py.iplot(fig)
La única conclusión definitiva que podemos extraer aquí es que el modelo trató muchos de los puntos de datos como valores atípicos.
La segunda función Prophet.plot_components podría ser mucho más útil en nuestro caso. Nos permite observar diferentes componentes del modelo por separado: tendencia, estacionalidad anual y semanal. Además, si proporcionamos información sobre vacaciones y eventos a nuestro modelo, también se mostrarán en este gráfico. Recalcar que este un caso muy especifico y solo tenemos un mes de historico(Mayo 2019)
m.plot_components(forecast);
La gráfica de la estacionalidad semanal lleva a la conclusión de que generalmente hay menos validaciones los sábados y domingos
Evaluación de calidad del pronóstico¶
Enfoque distinto a la primera parte
Para medir la calidad de nuestro pronóstico, necesitamos dividir nuestro conjunto de datos en la parte histórica , que es la primera y mayor porción de nuestros datos, y la parte de predicción , que se ubicará al final de la línea de tiempo. Eliminaremos los ultimos 10 días del conjunto de datos para usarlo más tarde como objetivo de predicción:
# Estos es como un train(493) y test(10)
prediction_size = 10
train_df = df89[:-prediction_size] # A excepcion de los ultimos 10 días
train_df.tail(n=5)
m = Prophet()
m.fit(train_df);
future = m.make_future_dataframe(periods=prediction_size)
future.tail(n=5)
Predecimos valores Prophet pasando las fechas para las que queremos crear un pronóstico. Si también proporcionamos las fechas históricas (como en nuestro caso), además de la predicción obtendremos un ajuste en la muestra para el historial. Llamemos al predict método del modelo con nuestro futuremarco de datos como entrada:
forecast = m.predict(future)
forecast.tail(n=5)
from fbprophet.plot import plot_plotly
import plotly.offline as py
import plotly.graph_objs as go
py.init_notebook_mode()
#m.plot(forecast);
fig = plot_plotly(m, forecast)
py.iplot(fig)
print(', '.join(forecast.columns))
Podemos ver que este marco de datos contiene toda la información que necesitamos, excepto los valores históricos. Necesitamos unir el forecast objeto con los valores reales ydel conjunto de datos original df. Para esto definiremos una función auxiliar que reutilizaremos más adelante:
def make_comparison_dataframe(historical, forecast):
"""Únimos a la historia con el pronóstico
El resultado contiene las columnas 'yhat', 'yhat_lower', 'yhat_upper' and 'y'.
"""
return forecast.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical.set_index('ds'))
# Habria que hacer un tratamiento mas exahustivo de datos, pero aqui se pretende explicar como funciona Prophet
cmp_df = make_comparison_dataframe(df89, forecast)
cmp_df.tail(n=5)
También vamos a definir una función auxiliar que usaremos para medir la calidad de nuestro pronóstico con medidas de error MAPE y MAE:
def calculate_forecast_errors(df89, prediction_size):
"""Calculate MAPE and MAE of the forecast.
Args:
df89: joined dataset with 'y' and 'yhat' columns.
prediction_size: number of days at the end to predict.
"""
# Make a copy
df89 = df89.copy()
df89['e'] = df89['y'] - df89['yhat']
df89['p'] = 100 * df89['e'] / df89['y']
# Recuerda que mantuvimos los valores de los últimos días de `prediction_size`
# para predecirlos y medir la calidad del modelo..
# Ahora recorta la parte de los datos para la que hicimos nuestra predicción. .
predicted_part = df89[-prediction_size:]
# Definimoa la función que promedia los valores de error absoluto sobre la parte predicha
error_mean = lambda error_name: np.mean(np.abs(predicted_part[error_name]))
# Ahora podemos calcular MAPE y MAE y devolver el diccionario de errores resultante.
return {'MAPE': error_mean('p'), 'MAE': error_mean('e')}
for err_name, err_value in calculate_forecast_errors(cmp_df, prediction_size).items():
print(err_name, err_value)
Como resultado, el error relativo de nuestro pronóstico (MAPE) es de aproximadamente 10%, y en promedio nuestro modelo está equivocado en 1638 validaciones (MAE).
from fbprophet.plot import plot_plotly
import plotly.offline as py
import plotly.graph_objs as go
py.init_notebook_mode()
fig = plot_plotly(m, forecast)
py.iplot(fig)
Visualizaciónes¶
Creemos nuestra propia visualización del modelo construido por Prophet. Comprenderá los valores reales, el pronóstico y los intervalos de confianza.
Primero, trazaremos los datos por un período de tiempo más corto para que los puntos de datos sean más fáciles de distinguir. En segundo lugar, mostraremos el rendimiento del modelo solo durante el período que predijimos, es decir, los últimos 10 días. Parece que estas dos medidas deberían darnos una trama más legible.
def show_forecast(cmp_df, num_predictions, num_values, title):
"""Visualize the forecast."""
def create_go(name, column, num, **kwargs):
points = cmp_df.tail(num)
args = dict(name=name, x=points.index, y=points[column], mode='lines')
args.update(kwargs)
return go.Scatter(**args)
lower_bound = create_go('Lower Bound', 'yhat_lower', num_predictions,
line=dict(width=0),
marker=dict(color="gray"))
upper_bound = create_go('Upper Bound', 'yhat_upper', num_predictions,
line=dict(width=0),
marker=dict(color="gray"),
fillcolor='rgba(68, 68, 68, 0.3)',
fill='tonexty')
forecast = create_go('Forecast', 'yhat', num_predictions,
line=dict(color='rgb(31, 119, 180)'))
actual = create_go('Actual', 'y', num_values,
marker=dict(color="red"))
# In this case the order of the series is important because of the filling
data = [lower_bound, upper_bound, forecast, actual]
layout = go.Layout(yaxis=dict(title='Validations'), title=title, showlegend = False)
fig = go.Figure(data=data, layout=layout)
py.iplot(fig, show_link=False)
show_forecast(cmp_df, prediction_size, 100, 'New passengers')
Tenemos un valor bajo de MAPE, puede ser debido a que ajusta demasiado bien es cierto que tenemos pocos datos.
Observamos como ejemplo en la muestra actual como el valor era de 15372 pero nos ha predicho que sera de 18k
df89[df89['ds']=='2019-05-27']
Podemos ajustar mas los datos con BOX-COX¶
import stats
from scipy import stats
def inverse_boxcox(y, lambda_):
return np.exp(y) if lambda_ == 0 else np.exp(np.log(lambda_ * y + 1) / lambda_)
train_df2 = train_df.copy().set_index('ds')
train_df2['y'], lambda_prophet = stats.boxcox(train_df2['y'])
train_df2.reset_index(inplace=True)
Creamos un nuevo Prophet modelo y repetimos el ciclo de ajuste de predicción que ya hemos hecho anteriormente:
m2 = Prophet()
m2.fit(train_df2)
future2 = m2.make_future_dataframe(periods=prediction_size)
forecast2 = m2.predict(future2)
for column in ['yhat', 'yhat_lower', 'yhat_upper']:
forecast2[column] = inverse_boxcox(forecast2[column], lambda_prophet)
cmp_df2 = make_comparison_dataframe(df89, forecast2)
for err_name, err_value in calculate_forecast_errors(cmp_df2, prediction_size).items():
print(err_name, err_value)
Entonces, definitivamente podemos indicar un minimo aumento en la calidad del modelo
Finalmente, tracemos nuestro rendimiento anterior con los últimos resultados de lado a lado. Tenemos en cuenta que usamos prediction_size para el tercer parámetro para acercar el intervalo que se predice:
show_forecast(cmp_df, prediction_size, 100, 'No transformations')
show_forecast(cmp_df2, prediction_size, 100, 'Box–Cox transformation')
Vemos que el pronóstico de los cambios semanales en el segundo gráfico está mucho más cerca de los valores reales ahora
- Conclusion:
Todavía depende del científico de datos explorar los resultados del pronóstico, ajustar los parámetros del modelo y transformar los datos cuando sea necesario.
Sin embargo, esta biblioteca es fácil de usar y fácilmente personalizable. La capacidad única de tener en cuenta los días anormales que el analista conoce de antemano puede marcar la diferencia en algunos casos.
En general, la biblioteca Prophet vale la pena que sea parte de nuestra caja de herramientas analíticas.